//EllipsoidModel（椭球模型）是用于模拟天体形状的数学模型，特别在描述地球形状时具有重要作用。以下是关于EllipsoidModel的一些详细解释：
//
//定义：
//EllipsoidModel可以视为一个定义了长短半轴的椭圆体，它用来模拟天体的形状，如太阳、行星和月亮等。在地球科学中，它常被用于模拟地球的形状。
//用途：
//与osg::CoordinateSystemNode（坐标系节点）一起使用，为场景中的对象关联一个坐标系统。
//实现经纬度与坐标之间的转换，这对于建立地球数据库和精确定位天体椭圆体模型的子节点非常有用。
//与WGS84的关系：
//在数字地球应用中，WGS84是一种常见的空间参考系统，它使用了一个特定的椭球模型来描述地球的形状。这个椭球模型具有特定的长半轴和扁率，使得它能够更准确地描述地球的实际形状。
//转换功能：
//EllipsoidModel能够实现经纬度（latitude和longitude）和高度（height）到世界坐标（X, Y, Z）的转换。这对于在GIS（地理信息系统）应用中进行数据的互通和融合至关重要，因为不同来源的地理数据可能采用不同的坐标系统。
//在OSG（OpenSceneGraph）中的应用：
//在OSG中，EllipsoidModel通过osg::WGS84Ellipsoid类实现，它提供了WGS84椭球模型的参数和方法。开发者可以使用这个类来进行坐标转换、计算法向量等操作。
//示例：
//在一些OSG示例中，如createEarth函数，会创建一个地球节点，并使用EllipsoidModel来定义地球的形状。然后，开发者可以使用这个模型来进行各种操作，如添加地理数据、计算位置等。
//总的来说，EllipsoidModel是一个强大的工具，它使得开发者能够更准确地模拟天体的形状，并在GIS应用中进行各种复杂的操作。
// Created by Administrator on 2024/6/10.
//

#ifndef COORDINATETRANSLATION_ELLIPSOIDMODEL_H
#define COORDINATETRANSLATION_ELLIPSOIDMODEL_H

#ifdef MATRIX_EIGEN
#include "eigen3/Eigen/Core"
#include "eigen3/Eigen/Geometry"
#else
#endif

namespace sgt {
#ifdef MATRIX_EIGEN
#define MATRIX4 Eigen::Matrix4d
#define Vector3 Eigen::Vector3d

#else
#define MATRIX osg::Matrixd
#define MATRIX osg::Matrixd

#endif
    const double PI   = 3.14159265358979323846;
    const double PI_2 = 1.57079632679489661923;
    const double PI_4 = 0.78539816339744830962;
    const double LN_2 = 0.69314718055994530942;
    const double INVLN_2 = 1.0 / LN_2;

    const float PIf   = 3.14159265358979323846f;
    const float PI_2f = 1.57079632679489661923f;
    const float PI_4f = 0.78539816339744830962f;
    const float LN_2f = 0.69314718055994530942f;
    const float INVLN_2f = 1.0f / LN_2f;


    const double WGS_84_RADIUS_EQUATOR = 6378137.0;
    const double WGS_84_RADIUS_POLAR = 6356752.3142;
    const double CGCS_2000_RADIUS_EQUATOR = 6378137.0;
    const double CGCS_2000_RADIUS_POLAR = 6356752.31414;


    inline float DegreesToRadians(float angle) { return angle*(float)PI/180.0f; }
    inline double DegreesToRadians(double angle) { return angle*PI/180.0; }

    inline float RadiansToDegrees(float angle) { return angle*180.0f/(float)PI; }
    inline double RadiansToDegrees(double angle) { return angle*180.0/PI; }
    ///它封装了用于模拟天文体（如太阳、行星、月球等）的椭球体。这个类提供了基本属性和方法，用于处理与椭球体相关的计算和转换
    class EllipsoidModel {
    public:
        ///默认使用WGS-84参数
        /// \param radiusEquator
        /// \param radiusPolar
        EllipsoidModel(double radiusEquator = WGS_84_RADIUS_EQUATOR, double radiusPolar = WGS_84_RADIUS_POLAR) :
                m_radiusEquator(radiusEquator), m_radiusPolar(radiusPolar) {
            computeCoefficients();
        }

        void setRadiusEquator(double radius) {
            m_radiusEquator = radius;
            computeCoefficients();
        }

        double getRadiusEquator() const { return m_radiusEquator; }

        void setRadiusPolar(double radius) {
            m_radiusPolar = radius;
            computeCoefficients();
        }

        double getRadiusPolar() const { return m_radiusPolar; }

        /// 将经纬高lat/long/hae 转为ECEF 坐标系下的XYZ
        /// \param latitude 纬度单位弧度
        /// \param longitude 经度单位弧度
        /// \param height 椭球高度，单位米
        /// \param X
        /// \param Y
        /// \param Z
        inline void convertLatLongHeightToXYZ(double latitude, double longitude, double height,
                                              double &X, double &Y, double &Z) const;

        /// 将ECEF坐标系下的位置转为经纬高lat/long/hae
        /// \param X
        /// \param Y
        /// \param Z
        /// \param latitude
        /// \param longitude
        /// \param height
        inline void convertXYZToLatLongHeight(double X, double Y, double Z,
                                              double &latitude, double &longitude, double &height) const;

        /// 获取 ENU站心坐标系 坐标系下的点转换到ECEF坐标系下的转换矩阵，该矩阵包含了旋转和平移
        /// \param latitude    站心坐标系原点所在的纬度，单位度
        /// \param longitude    站心坐标系原点所在的经度，单位度
        /// \param height       椭球高
        /// \param localToWorld 输出的旋转矩阵
        /// \attention 由ECEF坐标系到 ENU坐标系旋转矩阵是该矩阵的逆矩阵
        inline void computeLocalToWorldTransformFromLatLongHeight(double latitude, double longitude, double height,
                                                                  MATRIX4 &localToWorld) const;

        /// 获取 ENU站心坐标系 坐标系下的点转换到ECEF坐标系下的转换矩阵，该矩阵包含了旋转和平移
        /// \param X 站心坐标系原点所在的ECEF 向位置
        /// \param Y
        /// \param Z
        /// \param localToWorld 输出的旋转矩阵
        inline void computeLocalToWorldTransformFromXYZ(double X, double Y, double Z, MATRIX4 &localToWorld) const;

        /// 根据经纬度获取球面上一点为原点的的站心坐标系下的点转为ECEF点的旋转矩阵，不包括平移
        /// \param latitude     站心坐标系原点所在的纬度，单位度
        /// \param longitude    站心坐标系原点所在的经度，单位度
        /// \param localToWorld 输出的旋转矩阵
        inline void computeCoordinateFrame(double latitude, double longitude, MATRIX4 &localToWorld) const;

        inline Vector3 computeLocalUpVector(double X, double Y, double Z) const;

        // Convenience method for determining if EllipsoidModel is a stock WGS84 ellipsoid
        inline bool isWGS84() const {
            return (m_radiusEquator == WGS_84_RADIUS_EQUATOR && m_radiusPolar == WGS_84_RADIUS_POLAR);
        }

        // Compares two EllipsoidModel by comparing critical internal values.
        // Ignores _eccentricitySquared since it's just a cached value derived from
        // the _radiusEquator and _radiusPolar members.
        friend bool operator==(const EllipsoidModel &e1, const EllipsoidModel &e2) {
            return (e1.m_radiusEquator == e2.m_radiusEquator && e1.m_radiusPolar == e2.m_radiusPolar);
        }

    private:
        double m_radiusEquator;         //长半轴
        double m_radiusPolar;           //短半轴
        double m_eccentricitySquared;

    protected:
        /**
         * 据给定的长半轴和扁率计算椭球体的短半轴和其他系数
         */
        void computeCoefficients() {
            double flattening = (m_radiusEquator - m_radiusPolar) / m_radiusPolar;  //扁率（ƒ）
            m_eccentricitySquared = 2 * flattening - flattening * flattening;

//            // 计算短半轴
//            this.semiMinorAxis = semiMajorAxis * (1 - flattening);
//
//            // 额外计算椭球体的表面积和体积系数（如果需要）
//            // 椭球表面积系数（使用较准确的公式）
//            double surfaceAreaCoefficient = 4 * Math.PI * Math.pow(semiMajorAxis * semiMinorAxis * semiAxisOfRotation, 2.0 / 3.0);
//            // 这里假设semiAxisOfRotation为第三个轴的长度，但在天文体中可能并不适用
//            // 在实际应用中，椭球体可能没有明确的第三个轴长
//
//            // 椭球体积系数
//            double volumeCoefficient = (4.0 / 3.0) * Math.PI * semiMajorAxis * semiMinorAxis * semiAxisOfRotation;
//            // 注意：这里为了示例使用了semiAxisOfRotation，但在真正的椭球体模型中，可能只有两个轴（长半轴和短半轴）

        }
    };


} // sgt

#endif //COORDINATETRANSLATION_ELLIPSOIDMODEL_H
